function [y, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T,GAM0s,GAM1s,Cs,PSIs,PPIs,y0,eps_a, eps_u,complex_tol)

% NOTATION IS CONSISTENT WITH CAGLIARINI AND KULISH
% Solving Linear Rational Expectations Models with Predictable
% Structural Changes
%
% [y, S0, S1, S2, Sy, Sf, Sz, A, B] =
% LREAnt(n1,n2,s,T,GAM0s,GAM1s,Cs,PSIs,PPIs,y0,eps_a, eps_u,complex_tol)
% Takes the sequence of parameter matrices (GAM0s,GAM1s,Cs,PSIs,PPIs), the
% sequence of anticipated shocks (eps_a), the unanticipated shocks in the first period
% and the initial state (y0) to compute the state vector from period 1 to
% period T and the reduced-form matrices from time T+1 onwards
%
% Parameter matrices are stacked for periods 1 to T and then include the
% final system as the last set of parameters
%
% DIMENSIONS AND TIME PERIODS
% n1 is the number of non-expectational variables
% n2 is the number of expectational variables
% s is the number of horizons over which expectations appear in equations
% T is the number of time periods over which anticipated changes occurs
%
% TAKE SEQUENCE OF PARAMETERS
% GAM0s is n*(T+1) X n
% GAM1s is n*(T+1) X n
% Cs is n*(T+1) X 1
% PSIs is n*(T+1) X l  where l is the number of shocks
% PPIs is n*(T+1) X k  where k = n2*s
%
% INITIAL STATE VECTOR AND ANTICIPATED ADDITIVE SHOCKS AND UNANTICIPATED
% SHOCKS
% y0 is n X 1
% eps_a is l X T
% eps_u is l X 1

% THIS PROGRAM WAS WRITTEN FOR MATLAB VERSION 7.11 BY
%
%   ADAM CAGLIARINI AND MARIANO KULISH
%   ECONOMIC RESEARCH
%   RESERVE BANK OF AUSTRALIA
%   65 MARTIN PLACE
%   SYDNEY 2000, AUSTRALIA
%   adam.cagliarini@gmail.com
%   mkulish@gmail.com
%
%
% NOTE: Requires SMATS.M
% LAST MODIFIED: 3-June-2011
%%
k = s*n2;         % the dimension of expectational revisions
n = n1 + n2 + k;  % the dimension of the state variable
l = size(PSIs,2);  % the number of shocks

%% Check dimension assignments
 if k ~= size(PPIs,2);
     disp('The number of expectational revisions specified is different from the number of conditional expectations')
     disp('Check s, n2 or PPI')
     return
 end;
 
 if or(l~=size(eps_a,1),l~=size(eps_u,1))
     disp('The specified dimension of anticipated or unanticipated errors does not match the dimension implied by PSI')
     return
 end;
 
 if n~=size(GAM0s,2);
     disp('The specified dimension of the state vector does not match those implied by the parameter matrices')
     return
 end;
 
 if n~=size(y0,1);
     disp('Check dimension: y0 is not n x 1')
     return
 end;

 if l~=size(eps_u,1);
     disp('Check dimension: eps_u is not l x 1')
     return
 end;
 
%% Obtain the first and last set of parameter matrices
% The first set of parameters
GAM0o = GAM0s(1:n,:);
GAM1o = GAM1s(1:n,:);
Co = Cs(1:n);
PSIo = PSIs(1:n,:);
PPIo = PPIs(1:n,:);

% The last set of parameters
GAM0f = GAM0s(n*T+1:n*(T+1),:);
GAM1f = GAM1s(n*T+1:n*(T+1),:);
Cf = Cs(n*T+1:n*(T+1));
PSIf = PSIs(n*T+1:n*(T+1),:);
PPIf = PPIs(n*T+1:n*(T+1),:);

%% Remove the expectational errors from the first set of parameter matrices
[ind ignore] = find(sum(PPIo,2)==0);
GAM0ot = GAM0o(ind,:);
GAM1ot = GAM1o(ind,:);
Cot = Co(ind,:);
PSIot = PSIo(ind,:);

%% Compute the reduced-form matrices and the terminal condition w_{2,T}

% Compute the QZ decomposition of GAM0f and GAM1f and order so most
% explosive generalised eigenvalues are ordered last
[LAMBDA, OMEGA, Q, Z] = qz(GAM0f,GAM1f);
[LAMBDAS,OMEGAS,QS,ZS] = ordqz(LAMBDA,OMEGA,Q,Z,'udo');

lam = abs(diag(OMEGAS))./ abs(diag(LAMBDAS));

ind = find(lam>1);

nunstab = size(ind,1);

Q1 = QS(1:n-nunstab,:);
Z1 = ZS(:,1:n-nunstab);

Z2 = ZS(:,n-nunstab+1:n)';
Q2 = QS(n-nunstab+1:n,:);

m = rank(Z2);

Lambda11 = LAMBDAS(1:n-nunstab,1:n-nunstab);
Lambda12 = LAMBDAS(1:n-nunstab,n-nunstab+1:n);
Lambda22 = LAMBDAS(n-nunstab+1:n,n-nunstab+1:n);

Omega11 = OMEGAS(1:n-nunstab,1:n-nunstab);
Omega12 = OMEGAS(1:n-nunstab,n-nunstab+1:n);
Omega22 = OMEGAS(n-nunstab+1:n,n-nunstab+1:n);

M = Omega22\Lambda22;

W2T = (Lambda22 - Omega22)\(Q2*Cf);

if m==k;
    uniqueuptoT=1;
else;
    disp('Too many or too few terminal conditions');
    return;
end;


%% Sets up the A and B matrices and solves for the path;

A = [GAM0ot,zeros(n1+n2,(T-1)*n)];
B = Cot + GAM1ot*y0 + PSIot*(eps_a(:,1) + eps_u);

for i=2:T;
    A = [A; zeros(n,(i-2)*n), -GAM1s(n*(i-1)+1:n*i,:),GAM0s(n*(i-1)+1:n*i,:),zeros(n,(T-i)*n)];
    B = [B;Cs(n*(i-1)+1:n*i) + PSIs(n*(i-1)+1:n*i,:)*eps_a(:,i)];
end;

A = [A;zeros(m,(T-1)*n),Z2];
B = [B;W2T];

x = A\B;

y = zeros(n,T);

for i=1:T;
    y(:,i) = x((i-1)*n+1:i*n);
end;

if abs(sum(imag(x)))>complex_tol
    disp('Error: complex valued path')
    return
end
y = real(y); % Chops numerically insignificant remaining complex value parts.


%% Solve for the reduced-form matrices for the final system

[S0, S1, S2, S3, Sy, Sf, Sz, uniquea] = Smats(Cf, GAM0f, GAM1f, PSIf, PPIf);
